classdef GetMagnetometerBandwidth < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetMagnetometerBandwidth(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetMagnetometerBandwidth');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('nPoint', 50, @(x) x>0);
            obj.paramParser.addParameter('average', 100, @(x) x>0);
            obj.paramParser.addParameter('signal_amplitude', 0.01, @(x) x>0 && x < 0.1);
            obj.paramParser.addParameter('start_frequency', 10.0, @(x) x>0);
            obj.paramParser.addParameter('finish_frequency', 300.0, @(x) x>0);
            obj.paramParser.addParameter('order', 4, @(x) ismember(x, 1:8));
            obj.paramParser.addParameter('bandwidth', 10.19, @(x) x>0);
            obj.paramParser.addParameter('liaIdx', 3, @(x) ismember(x, 1:4));
            obj.paramParser.addParameter('f129', 250.0, @(x) isnumeric(x) && x>0);
            obj.paramParser.addParameter('f131', 75.0, @(x) isnumeric(x) && x>0);
            obj.paramParser.addParameter('fRange', 40, @(x) isnumeric(x) && x>0);
            obj.paramParser.addParameter('description', 'normal', @ischar);
        end
        
        %%
        function obj = preprocess(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;

            nmrg.UnlockTransverseField('x'); pause(0.5);
            nmrg.bx.setDriver( nmrg.ziSlave.LockInChannel(1).setOscillator(param.liaIdx) );
            nmrg.bx.setAC('frequency', 10, 'amplitude', param.signal_amplitude, 'offset', nmrg.SystemInfo.WorkingPoint.biasX);
        end
        
        %%
        function obj = dataAcquisition(obj)
            import MeasurementComponent.Signal.*
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            
            measurement_data.start = datetime('now');
            liaChannel = nmrg.ziSlave.LockInChannel(param.liaIdx).setup('order', param.order, 'bandwidth', param.bandwidth, 'phase', 0.0,  'isSinc', true);
            bxSig = Signal('bxAmplitude', liaChannel, 'reference', nmrg.bx); pause(5);
            
            bxSig.subscribe_signal(ziDemodTime.R_avg, [2, 1, 1]);
            bxSig.subscribe_signal(ziDemodTime.Theta_avg, [2, 1, 2]);
            measurement_data.sweeperData = bxSig.getParameterDependence(nmrg.ziSlave, ziSweepNode.Osc3Frequency, param.start_frequency, param.finish_frequency, ...
                                                                        'npoint', param.nPoint, 'averaging_sample', param.average,'isLog', false, 'sincfilter', 1);
            measurement_data.finish = datetime('now');
            obj.dataPackage.measurement_data = measurement_data;
        end
        
        %%
        function obj = analyze(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measurementData = obj.dataPackage.measurement_data;

            amp = param.signal_amplitude;
            f129 = param.f129;
            f131 = param.f131;
            fRange = param.fRange;
            swprData = measurementData.sweeperData;
            
            %-----------------------------------
            maxAmp = nmrg.SystemInfo.WorkingPoint.slopeX*amp*nmrg.SystemInfo.WorkingPoint.exportScaleAC/sqrt(2);
            data1 = swprData.getNodeData(1); fList = data1.x; dataR = data1.y;
            [xData, yData] = prepareCurveData(fList, dataR);
            
            excludeIndex = find( (f131-fRange<xData & xData<f131+fRange) | (f129-fRange<xData & xData<f129+fRange));
            excludedPoints = excludedata( xData, yData, 'Indices', excludeIndex);
            ft = fittype('a*g/sqrt(x^2+g*g)+c', 'independent', 'x', 'dependent', 'y' );
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.StartPoint = [maxAmp, min(yData), 20];
            opts.Exclude = excludedPoints;
            
            [fittingR.result, fittingR.gof] = fit( xData, yData, ft, opts );
            fittingR.dcValue = maxAmp;
            fittingR.amplitude = fittingR.result.a;
            fittingR.background = fittingR.result.c;
            fittingR.bandwidth = fittingR.result.g;
            fittingR.exclude = [f129, f131];
            fittingR.excludeRange = fRange;
            fittingR.gain129 = fittingR.result(f129)/amp;
            fittingR.gain131 = fittingR.result(f131)/amp;
            
            %-----------------------------------
            data2 = swprData.getNodeData(2); fList=data2.x; dataTheta=data2.y;
            [xData, yData] = prepareCurveData(fList, dataTheta);
            excludedPoints = excludedata( xData, yData, 'Indices', excludeIndex);
            ft = fittype( 'smoothingspline' );
            opts = fitoptions( 'Method', 'SmoothingSpline' );
            opts.Exclude = excludedPoints;
            
            [fittingTheta.result, fittingTheta.gof] = fit( xData, yData, ft, opts );
            fittingTheta.exclude = [f129, f131];
            fittingTheta.excludeRange = fRange;
            fittingTheta.phase129 =fittingTheta.result(f129)/pi*180;
            fittingTheta.phase131 =fittingTheta.result(f131)/pi*180;
            
            %-----------------------------------
            analysis_result.R = fittingR;
            analysis_result.Theta = fittingTheta;
            obj.dataPackage.analysis_result = analysis_result;
        end
        
        %%
        function obj = postprocess(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;
            
            %-----------------------------------
            nmrg.bx.disable();
            nmrg.ziSlave.LockInChannel(param.liaIdx).setup('isSinc', false);
            nmrg.LockTransverseField('x');
            
            %-----------------------------------
            summary.start = datestr(measuredData.start);
            summary.polarization = 'x';
            summary.amplitude = param.signal_amplitude;
            summary.start_frequency = param.start_frequency;
            summary.finish_frequency = param.finish_frequency;
            summary.filter_order = param.order;
            summary.filter_bandwidth = param.bandwidth;
            summary.filter_sinc = 'on';
            summary.gain129 = analysis.R.gain129;
            summary.gain131 = analysis.R.gain131;
            summary.phase129 = analysis.Theta.phase129;
            summary.phase131 = analysis.Theta.phase131;
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            nmrg.SystemInfo.MagnetometerBandwidth = summary;
            
            %-----------------------------------
            getSession;
            sess.addData('GetMagnetometerBandwidth', obj.dataPackage, param.description);
        end
        
        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            swprData = measuredData.sweeperData;
            
            if nargin == 1
                fig = figure('Name', sprintf('Bandwidth: [%s]', datestr(measuredData.start)), 'WindowStyle', 'Docked'); 
                ax1=subplot(2, 1, 1, 'Parent', fig);
                ax2=subplot(2, 1, 2, 'Parent', fig);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);
                ax1=subplot(2, 1, 1,  'Parent', viewer.gui.uiHandle.panelPlot);
                ax2=subplot(2, 1, 2,  'Parent', viewer.gui.uiHandle.panelPlot);
            end
            
            data1 = swprData.getNodeData(1); fList=data1.x; amplitude=data1.y;
            data2 = swprData.getNodeData(2); phase = data2.y;
            fList0 = linspace(0, max(fList), 101);
            
            plot(ax1, fList, amplitude, 'b.', fList0, analysis.R.result(fList0), 'r-');
            xlim(ax1, [0 max(fList)]); grid(ax1, 'on');
            xlabel(ax1, 'Driving Frequency (Hz)'); ylabel(ax1, 'Magnetometer Amplitude (V)');
            title(ax1, sprintf('Signal Amplitude %3.2f V, BW = %3.1f Hz, Amp. = %5.3fV, Bg. = %5.3fV', ...
                                summary.amplitude, analysis.R.bandwidth, analysis.R.amplitude, analysis.R.background));

            plot(ax2, fList, unwrap(phase)/pi*180, 'b.', fList0, analysis.Theta.result(fList0)/pi*180, 'r-');
            xlim(ax2, [0 max(fList)]); grid(ax2, 'on');
            xlabel(ax2, 'Driving Frequency (Hz)'); ylabel(ax2, 'Magnetometer Phase (deg)');
            title(ax2, sprintf('g_{129}=%3.1f, \\theta_{129} = %3.2f^{\\circ}; g_{131}=%3.1f, \\theta_{131} = %3.2f^{\\circ}', ...
                analysis.R.gain129, analysis.Theta.phase129, analysis.R.gain131, analysis.Theta.phase131));
        end
    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                startTime = arrayfun(@(t) t.start, res.summary, 'UniformOutput', false);
                if isempty(p.Results.viewer)
                    f = figure('Name', sprintf('Bandwidth'));
                    ax1 = subplot(2, 1, 1, 'Parent', f);
                    ax2 = subplot(2, 1, 2, 'Parent', f);
                else
                    ax1 = subplot(2, 1, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax2 = subplot(2, 1, 2, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                
                nlen = length(res.measuredData);
                rData = arrayfun(@(d) d.sweeperData.getNodeData(1), res.measuredData);
                phiData = arrayfun(@(d) d.sweeperData.getNodeData(2), res.measuredData);
                
                plot(ax1, reshape([rData.x], [], nlen), reshape([rData.y], [], nlen), '.-'); grid(ax1, 'on');
                xlabel(ax1, 'Driving Frequency (Hz)'); ylabel(ax1, 'Magnetometer Amplitude (V)');
                legend(ax1, startTime, 'Location', 'bestoutside', 'Interpreter', 'tex');
                
                plot(ax2, reshape([phiData.x], [], nlen), unwrap(reshape([phiData.y], [], nlen))/pi*180, '.-'); grid(ax2, 'on');
                xlabel(ax2, 'Driving Frequency (Hz)'); ylabel(ax2, 'Magnetometer Phase (deg)');
                legend(ax2, startTime, 'Location', 'bestoutside', 'Interpreter', 'tex');
            end
        end
    end
end

